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METHOD FOR SECOND-ORDER ELLIPTIC EQUATIONS 
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Abstract. The weak Galerkin finite element method is a novel numerical method that was first 
proposed and analyzed by Wang and Ye in |23l for general second order elliptic problems on triangular 
meshes. The goal of this paper is to conduct a computational investigation for the weak Galerkin 
method for various model problems with more general finite element partitions. The numerical results 
confirm the theory established in |23j . The results also indicate that the weak Galerkin method is 
efficient, robust, and reliable in scientific computing. 
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1. Introduction. In this paper, we are concerned with computation and numer- 
ical accuracy issues for the weak Galerkin method that was recently introduced in j23| 
for second order elliptic equations. The weak Galerkin method is an extension of the 
standard Galerkin finite element method where classical derivatives were substituted 
by weakly defined derivatives on functions with discontinuity. The weak Galerkin 
method is also related to the standard mixed finite element method in that the two 
methods are identical for simple model problems (such as the Poisson problem). But 
they have fundamental differences for general second order elliptic equations. The 
goal of this paper is to numerically demonstrate the efficiency and accuracy of the 
weak Galerkin method in scientific computing. In addition, we shall extend the weak 
Galerkin method of J23' from triangular and tetrahedral elements to rectangular and 
cubic elements. 

For simplicity, we take the linear second order elliptic equation as our model 
problem. More precisely, let ft be an open bounded domain in R'', d = 2,3 with 
Lipschitz continuous boundary dil. The model problem seeks an unknown function 
u = u{x) satisfying 

-V ■ (AVu) + (3 -Vu + ju = / inn, 
^ ' ' u = g on on, 

where A G [L°^ (n)]'^'"^ , (3 G [L°°{n)]'^, and 7 G L°°{n) are vector- and scalar-valued 
functions, as appropriate. Furthermore, assume that ^ is a symmetric and uniformly 
positive definite matrix and the problem (jl.ip has one and only one weak solution in 
the usual Sobolev space H^{rt) consisting of square integrable derivatives up to order 
one. / and g are given functions that ensure the desired solvability of 

Throughout the paper, we use || • || to denote the standard norm over the 
domain fl, and use bold face Latin characters to denote vectors or vector- valued 
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functions. The paper is organized as follows. In Section 2, the weak Galerkin method 
is introduced and an abstract theory is given. In particular, we prove that certain 
rectangular elements satisfy the assumptions in the abstract theory, and thus establish 
a well-posedness and error estimate for the corresponding weak Galerkin method with 
rectangular meshes. In Section 3, we present some implementation details for the weak 
Galerkin elements. Finally in Section 4, we report some numerical results for various 
test problems. The numerical experiments not only confirm the theoretical predictions 
as given in the original paper |23) , but also reveal new results that have not yet been 
theoretically proved. 

2. The Weak Galerkin Method. Let 7/i be a shape-regular, quasi-uniform 
mesh of the domain f2, with characteristic mesh size h. In two-dimension, we consider 
triangular and rectangular meshes, and in three-dimension, we mainly consider tetra- 
hedral and hexahedral meshes. For each element K € Th, denote by Kg and dK the 
interior and the boundary of respectively. Here, K can be a triangle, a rectangle, 
a tetrahedron or a hexahedron. The boundary dK consists of several "sides" , which 
are edges in two-dimension or faces(polygons) in three-dimension. Denote by J-^ the 
collection of all edges/faces in Th- 

On each if G 7h, let Pj{Ko) be the set of polynomials on Kq with degree less 
than or equal to j, and Qj{K) be the set of polynomials on Kq with degree of each 
variable less than or equal to j. Likewise, on each F E Fy^ Pi{F) and Qi{F) are 
defined analogously. Now, define a weak discrete space on mesh Th by 



Observe that the definition of Sh does not require any form of continuity across 
element or edge/face interfaces. A function in Sh is characterized by its value on the 
interior of each element plus its value on the edges/ faces. Therefore, it is convenient 
to represent functions in Sh with two components, v — {wq, ^^b}, where denotes the 
value of V on all KqS and Vb denotes the value of v on Fh- 

We further define an L? projection from H^{il) onto Sh by setting QhV = 
{Qov, Qbv}, where Qov\k is the local projection of v in Pj{Ko), for K e Th, 
and Qbv\F is the local projection in Pi{F), for F G Fh- 

The idea of the weak Galerkin method is to seek an approximate solution to 
Equation in the weak discrete space Sh- To this end, we need to introduce a 

discrete gradient operator on Sh- Indeed, this will be done locally on each element 
K. Let Vr{K) be a space of polynomials on K such that [Pr{K)]'^ C Vr{K); details 
of Vr{K) will be given later. Let 



A discrete gradient of Vh = {vo,Vb} G Sh is defined to be a function V^f/i G S/i such 
that on each K G Th, 



where n is the unit outward normal on dK- Clearly, such a discrete gradient is always 
well-defined. 

Denote by (•,•) the standard L^-inner product on il. Let S'fl be a subset of Sh 
consisting of functions with vanishing boundary values. Now we can write the weak 



Sh ^{v ■- v\ko G Pj{Ko) or Qj{Ko) for all K G Th, 
v\f G Pi{F) or Qi{F) for all F G Fh}- 



S/. = {q e [L''in)f ■- q\K G Vr{K) for aU K G %}- 



(2.1) 




for all q G Vr{K), 
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Galerkin formulation for Equation (|l.ip as follows: find Uh = {wqi ui,} G Sh such that 
Ub = QbU on each edge/face F C dil and 

(2.2) (^VdUft, VdV,,) + (/3 ■ VdU,,, Wo) + (7M0, Wo) = (/, «o) 

for all Vfi — {vQ,Vb} G S'^^. For simplicity of notation, we introduce the following 
bilinear form 

(2.3) a{uh, Vh) = [AV dUh,^ dVh) + (^ • ^ dUh, vq) + {'juo, Vq). 

The spaces Sh and S/i can not be chosen arbitrarily. There are certain criteria 
they need to follow, in order to guarantee that Equation p.2p provides a good ap- 
proximation to the solution of Equation (jl.ip . For example, Uh has to be rich enough 
to prevent from the loss of information in the process of taking discrete gradients, 
while it should remain to be sufficiently small for its computational cost. Hence, we 
would like to impose the following conditions upon Sh and S/j: 

(PI) For any Vh € Sh and K G Th, V^w/ilif = if and only if vq = Vb = constant 
on K. 

(P2) For any w G iJ™+^(r2), where < m < j + 1, we have 

\\^d{Qhw)~Vw\\<Ch'^\\w\\^+i, 

where and in what follows of this paper, C denotes a generic constant inde- 
pendent of the mesh size h. 

Under the above two assumptions, it has been proved in [23] that Equation (|2.2|) 
has a unique solution as long as the mesh size h is moderately small and the dual of 
(|l.ip has an iJ^+^-regularity with some s > 0. Furthermore, one has the following 
error estimate: 

W^diuh - Qhu)\\ < C {h'+-'\\f - QofW + , 
\\uo~Qou\\<C{h'+'\\f-Qof\\+h"'+'\\u\\r„+,), 

for any < m < j + 1, and s > is the largest number such that the dual of Equation 
(jl.ip has an iJ^+'^-regularity. 

There are several possible combinations of Sh and E/i that satisfy Assumptions 
(PI) and (P2). Two examples of triangular elements have been given in [53], which 
are 

1. Triangular element {Pj{Ko), Pj{F), RTj{K)) for j > 0. That is, in the def- 
inition of Sh, we set I ~ j. And in the definition of E/j, we set r — j and 
VriK) to be the jth order Raviart-Thomas element RTj{K) fST. 

2. Triangular element {Pj{Ko), P]+i{F), {Pj+i{K))'^) for j > 0. That is, in 
the definition of Sh, we set I = j + I. And in the definition of S/i, we set 
r = j + 1 and Vr{K) = {Pj+i{K)Y, or in other words, the {j + l)st order 
Brezzi-Douglas-Marini element [^ . 

The rest of this section shall extend this result to rectangular elements. An ex- 
tension to three-dimensional tetrahedral and hexahedral elements is straightforward. 

2.1. Weak Galerkin on Rectangular Meshes. Consider the following two 
type of rectangular elements: 
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1. Rectangular element {Qj{Ko), Qj{F), RTj{K)) for j > 0. That is, in the 
definition of S*^, we set I — j. And in the definition of E/j, we set r = j and 
Vr{K) to be the jth order Raviart-Thomas element RTj{K) on rectangle K. 

2. Rectangular element {Pj{Ko), Pj+i{F), BDMj+i{K)) for j > 0. That is, in 
the definition of Sh, we set / — j + 1. And in the definition of E/i, we set 
r = j + 1 and Vr{K) to be the {j + l)st order Brezzi-Douglas-Marini element 
BDMj+i{K) on rectangle K. 

Denote by Qij{K) the space of polynomials with degree in x and y less than or equal 

-d/dy 
d/dx 



to i and j, respectively, and curl 



. It is known that 



BDMj+i{K) = 



span 



{curl (a;-'^^y), curl (xy-''*'^)} , 



and 6:mY{RTj{K)) = 2(j + l)(j + 2), 6:in\{BDMj+i{K)) = (j + 2)(j + 3) + 2. The 
degrees of freedom for RTj{K) are: 

/ (q • n)w ds, for ah w € Qj{F), F e T D dK, 

J F 

q • pdcc, for all p G Qj-ij{K) x Qjj_i{K). 

The degrees of freedom for BDMj^i{K) are 

/ (q • n) w ds, for all w ^ Pj+i{F), F e n dK, 

JF 



q • Tpdx, 



for all p e 



It is also well-known that on each rectangle K G Th and each edge F E ThD dK, 



(2.5) 



\7 ■RTj{K) = Qj{Ko) 
V ■ BDM,+i{K) ^ Pj{Ko) 



RT,{K)-n\F^Q.j{F), 
BDM,+i{K)-n\F=Pj+i{F). 



Next, we show that the two set of elements defined as above satisfy Assumptions 
(PI) and (P2). 

Lemma 2.1. For the two type of rectangular elements given in this subsection, 
the Assumption PI holds true. 

Proof. If vq = Vh = constant on K, then clearly VdW;i|if vanishes since the 
right-hand side of (|2.1[) is zero from the divergence theorem. Now let us assume that 
Vd'yft.lif = 0. By (|2.ip and using integration by parts, we have for all q G RTj{K) or 
BDM,+,iK), 







(2.6) 



wqV ■ qdx 



K 



Wbq • n ds 



dK 



I (wfc - wo)q • nds + / {yvQ)-q^dx. 

JdK Jk 



We first consider the element {Qj{KQ), Qj{F), RTj{K)). If j = 0, then vq is a 
constant on Kq and clearly Vup = 0. If j > 0, take q such that /p(q • n)w ds — for 
all w G Qj{F) and let it traverse through all degrees of freedom defined by q-pdx, 
for p e Qj^i^j{K) X Qjj_i{K). Since {vi, - vo)\f e Qj{F) and Vuq £ Q]-i,j{K) x 
Qjj-i{K), Equation (j2.6p gives Vuq = 0, which implies that vq is a constant on Kq. 
Now Equation (|2.6p reduces into 

/ {vb - vo)q • n ds = 0, for ah q G RTj{K). 

JdK 

Next, since (vb - vo)\f e (9j(F) = RTj{K) ■ n\F for all F e n dK, by letting q 
traverse through all degrees of freedom on dK, we have Vb — vo ~ on all F. This 
implies Vb = vo = constant in _fC. 

For the {Pj{Ko), Pj+i(F), BDMj+i{K)) element, using the same argument as 
in the previous case, and noticing that Vwq G {Pj-i{K)y , {vb — wo)|_f & Pj+i{F) = 
BDMj+i{K)-n\p for all F S Fhr\dK , we can similarly prove that = wq = constant 
in iC. □ 

Lemma 2.2. For the two type of rectangular elements given in this subsection, 
the Assumption (P2) holds true. 

Proof. Let w G i?'"(r2), < m < j + 1. For any q e S?, and £ Th, by ([23]) 
and the definition of projections, we have 

' (y dQhw) ■ q_dx ^ - / ((5ow)(V • q) da; + / (9bw(q-n)(is 
= — w(V • q) da; + / w(q • n) ds 

JK JdK 

— I (Vw) • qda;. 

JK 

In other words, on each K e 7/t, ^dQhW is the L'^ projection of Vw onto RTj{K) or 
BDMj^i{K). Thus, the Assumption (P2) follows immediately from the approxima- 
tion properties of the projection, and the fact that both RTj{K) and BDMj^i{K) 
contains the entire polynomial space [Pj{K)]'^. □ 

Using Lemma [23] and Lemma [221 one can derive the error estimate ()2.4|) for the 
rectangular elements by following the argument presented in [23]. Details are left to 
interested readers as an exercise. 

3. Computation of Local Stiffness Matrices. Similar to the standard Galerkin 
finite element method, the weak Galerkin method (|2.2p can be implemented as a ma- 
trix problem where the matrix is given as the sum of local stiffness matrices on each el- 
ement K € Th- Thus, a key step in the computer implementation of the weak Galerkin 
is to compute element stiffness matrices. The goal of this section is to demonstrate 
ways of computing element stiffness matrices for various elements introduced in the 
previous sections. 

For a given element K G Th, let (/)o,i, i = l,...,No, be a set of basis func- 
tions for Pj{Ko) or Qj{Ka), and (pb.i, i — 1, . . . ,Nb, be a set of basis functions for 

^FedKnJ^h ■^'(-^) °^ J^FedKnJ^h Qi(^) ■ ^o*"^ ^^^^ i'f'b^t} is the union of basis func- 
tions from all edges/faces of element K. Then every Vh = {vo,Vb} £ Sh has the 
following representation in K: 



Vh\K 



On each K, the local stiffness matrix for Equation (|2.2p can thus be written as 
a block matrix 



(3.1) 



Mk^ 



'Mq.,0 Mo.b 
Mb.o Mb.b 



where Mq.o is an Nq x Nq matrix, Afo,6 is an A^o x matrix, Mb.o is an Nb x iVo 
matrix, and Mb.b is an Nb x Nb matrix. These matrices are defined, respectively, by 



Mofi = [a(0oj,0o,i)K] 
MbM = [a(0oj, 06,i)k] 



Mo,b = [ai4>b,j,(l}o,i)K] 
Mbfi = [a[cf)b,j,(j)b,i)K] 



where the bilinear form a(-, •) is defined as in (|2.3p . and i, j are the row and column 
indices, respectively. 

To compute each block of Mk, we first need to calculate the discrete gradient 
operator V^. For convenience, denote the local vector representation of Vh\K by 





Vo.l 




Vb,l 




Vq.2 




Vb.2 


Xo = 




y.b = 






yo.No_ 




yb,Nh_ 



Let Xii i = ■ ■ 1 ^v, be a set of basis functions for Vr{K). Then, for every q/i € T,h, 
its value on K can be expressed as 



Nv 



Similarly, we denote the local vector representation of (IhIk by 



q = 



Then, by the definition of the discrete gradient (|2.ip . given Vh\K, we can compute the 
vector form of V^w/i on K by 



(3.2) 



where the Ny x Ny matrix Dk, the Ny x A^o matrix Zk, and the Ny x Nb matrix 
Tk are defined, respectively, by 



(3.3) 



Dk = 



Ik Xi ■ Xi dx 
Jk Xnv ■ Xi dx 

Sk^"^ ■Xi)<pQ.idx 



Ik Xi ■ Xnv 
Ik Xnv ■ Xnv dx_ 

Ixi^ ■ Xi^o.No dx 
• Xnv)^o,No dx_ 



and 



Notice that Dk is a symmetric matrix. 

Once the matrices D^, Zx and are computed, we can use (j3.2l) to calculate 
the weak gradient of basis functions (jjQ i and on It is not hard to see that 



(3.4) 



No 



where e^" and e^' are the standard basis for the Euclidean spaces and M^'' 
respectively, such that its i-th entry is 1 and all other entries are 0. 
Define matrices 



(3.5) 



[(/3 • XjAo,i)K_ 



Bk 
Ck 



[(70o,j,'/'o,i) 



where (•, ■)k denote the standard inner-product on L'^{K) or [L'^{K)]'^, as appropriate. 
Clearly, Ak is an Ny x Ny matrix, Bk is an A'o x Ny matrix, and Ck is an A^o x Nq 
matrix. Then, an elementary matrix calculation shows that the local stiffness matrix 
Mk for Equation (|2.2p can be expressed in a way as specified in the following lemma. 



Lemma 3.1. The local stiffness matrix Mk defined in iS.l\) can be computed by 
using the following formula 



(3.6) 



^^0,0 — Z]^Dj^AkD^Zk — BkDj^Zk + 
Mo,6 = -Z]iD^'AKD^^TK + BkD^^Tk. 
Mb,o = -T'i,Dj,'AkDj^'Zk + Tl^D-^'B'x, 
Mb,b = T'j,D-'AkD]^^Tk, 



C. 



K, 



where the superscript t stands for the standard matrix transpose. 

For the Poisson equation — Au — /, we clearly have Ak — Dk and Bk — 0, 
Ck ~ 0. Since Dk is symmetric, the local stiffness matrix becomes 



(3.7) 



M 



K 



^k^kZk -Z'^Dj/Tk 
-^^kD^^Zk T^^D^^Tk ^ 



In the rest of this section, we shall demonstrate the computation of the element 
stiffness matrix AIk with two concrete examples. 

3.1. For the Triangular Element (Po(K), Po(F), RTo(K)). Let be a tri- 
angular element in Th- We consider the case when j = I = Q and Vr{K) being the 
lowest order Raviart-Thomas element. In other words, the discrete space Su con- 
sists of piecewise constants on the triangles, and piecewise constants on the edges 
of the mesh. In this case, the discrete gradient is defined by using the lowest order 
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Raviart-Thomas element on the triangle K. Clearly, we have A^o = 1) -^6 = 3 and 

Nv = 3. 

Let Vi = {xi,yi), i = 1,2,3, be the vertices of the triangle K and Cj be the edge 
opposite to the vertex Vi. Denote by \ei\ the length of edge and \K\ the area of the 
triangle K . We also denote by rij and t, the unit outward normal and unit tangen- 
tial vectors on Cj, respectively. Here should be in the positive (counterclockwise) 
orientation. If edge e, goes from vertex vj to Vk and K stays on the left when one 
travels from Vj to Vk, then it is not hard to see that 



U = 



'u.i 


1 




Xj 


ti,2 


"N 


Vk 


-yj. 



n,; 





1 


Vk - Vj 


ni,2 


"R 





3.1.1. Approach I. One may use the following set of basis functions for the 
weak discrete functions on K: 



(3.8) 

and 
(3.9) 



<t>o, 



on Cj 
otherwise 



for i = 1,2,3, 



Xi = 



2\K\ 



Xi 

Vi 



for i = 1, 2, 3. 



Notice that forms the standard basis for the lowest order Raviart-Thomas clement, 
for which the degrees of freedom are taken to be the normal component on edges. 
Indeed, Xi satisfies 



1 for i = j, 

for i 7^ j. 

It is straight forward to compute that, for the above defined basis functions. 



Zk = 



k'2| 



|ei| 







|e2| 








leal 



The computation of Dk is slightly more complicated, but it can still be done without 
much difficulty, especially with the help of symbolic computing tools provided in 
existing software packages such as Maple and Mathematica. For simplicity of notation, 
denote 



— l^il 
Uj = l^il + |ej| 

^123 = |eip + |e2| 
Then, it can be verified that 



for 1 < i < 3, 

for 1 < i, J < 3 and i ^ j, 



Dk = 



1 



4&\K\ 



(3.10) 



|eip(3;23- 
|ei||e2| (^12 
|ei||e3| (Zi3 



|e3r- 

■3/3) 
■3Z2) 



|ei||e2| (/12 
|e2|'(3^i3 
|e2||e3| (^23 



3^3) 

■^2) 
■3/1) 



|ei||e3|(Zi3-3/2)' 
|e2||e3| (/23 - 3Zi) 
|e3|' (3Zi2 - h) 



AS,\K 



■Tk 



3^23 — ^1 ^12 — 3/3 lis — 3Z2 

Z12 — 3Z3 3/13 — I2 I23 — 3Zi 

^13 — 3Z2 ^23 — 3^1 3Zi2 — h 



We point out that, the value of given as in (I3.10p agrees with the one presented 
in [3] . A verification of the formula (I3.10|) can be carried out by using the following 
fact 



' ' 2 



1 1 1 

Xi X2 

yi y2 ys 



{xj - xkf + {yj - Vkf 



1, 2, 3, j 7^ fc, j and k different from i. 



In computer implementation, it is convenient to use a form for the local matrix that 
can be expressed by using only edge lengths, as the one given by (|3.10p . 

In addition, using symbolic computing tools and the law of sines and cosines, we 



can write D as follows: 



D 



K 



T, 



K 



161X1 



'123 



21X1 



■ 2li 

h — ha 



h — h2 



h — ha 
h ~ hs 
2I3 



K 



Thus, to compute the local stiffness matrix it suffices to calculate Ak, Bk 
and C/f as given in p.Sp . and then apply Lemma l3.ll Notice that these three matrices 
depend on the coefficients A, /3 and 7, and quadrature rules may be employed in the 
calculation. However, for the simple case of the Poisson equation — Au = /, we see 
from (l3?7l) that 



144|A: 



Mi 



Ofc 



bO 



16|g| 

h 



23 



1 1 
1 1 
1 1 



21X1 





-iS\K\ 


-48\K\ 


-48\K\ 












2h 


h ^ h2 


h ~ hs 


^3 


-h2 


2I2 


h — 123 


h 


-hs 


h — I23 


2h 



3.1.2. Approach II. We would like to present another approach for computing 
the local stiffness matrix Mk in the rest of this subsection. Observe that a set of 
basis functions for the local space Vr{K) can be chosen as follows 



(3.11) 



Xi 



X2 



X3 



X — X 

y-y 



where {x = (xi +X2 +X3)/3,y = (yi + j/2 +2/3)73) is the coordinate of the barycenter 
of K. Note that both components of X3 have mean value zero on K. For the weak 
discrete function on K, we use the same set of basis functions as given in p.8p . It is 
not hard to see that 



Dk = |X| 



Hf- 

36 . 



Zk = 



and 



Tk = 



J/3 - 2/2 yi- y3 

X2 - X3 X3- Xi 

Ml Ml 

3 3 
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y2 



Xi 






2|X1 

yi 

X2 



2\K\ 
3 



Next, we use the formula (I3.5P to calculate the matrices A^, and Ck for the new 
basis (|3.1ip . Finally, we calculate the local stiffness matrix Mj^ by using the formula 
provided in Lemma [3. II 

Since the set of basis functions for the weak discrete space are the same in Ap- 
proaches I and II, the resulting local stiffness matrix AI^ would remain unchanged 
from Approaches I and II. The set of basis functions (I3.11|) is advantageous over the 
set (j3.9p in that the matrix Dk is a diagonal one whose inverse in trivial to compute. 



3.2. For the Cubic Element (Qo(K), Qo(F), RTo(K)). Let K = [0,a] x 

[0, 6] X [0, c] be a rectangular box where a, b, c are positive real numbers. We consider 
the three-dimensional cubic element, for which the discrete space Sh consists of piece- 
wise constants on Kq and piecewise constants on the faces of K. The space for the 
discrete gradient is the lowest order Raviart-Thomas element on K. We clearly have 
Nq = I, Nb = 6 and Ny = 6. 

Denote the six faces F^, i = 1, . . . , 6 by 



Fs : y = 0, 
: z = 0. 



F2 : X = a, 
Fa ■■ y^h, 
Fe : z^c. 



Note that the volume of K is given by \K\ — abc and the normal direction to each 
face is given by 





"-1" 




"1" 




" " 




'0 




" " 




"0' 


ni = 





, na = 





, ng = 


-1 


, n4 = 


1 


, ng = 





, ng = 



























-1 




1 



We adopt the following set of basis functions for the weak discrete space on K 

f 1 onF, 

00,1 = 1, 0b,i = < for I = 1,...,6, 

otherwise 



and 





a 




~ x' 

a 









"0" 









'0' 


Xi = 





, X2 = 





. X3 = 




, Xi 


y 

b 







. Xe = 



























^ - 1 

_ c 




z 
_ c_ 



Clearly, each Xi satisfies 



X,-nj\F, 



It is not hard to compute that 



\K\ 



2 

-1 







-1 
2 










-1 

2 





for i = j, 
for i ^ j. 



D 



K 



2 

W\ 



2 1 

1 2 

2 1 

1 2 

2 1 

1 2 



10 



and 



Z 



K 



'be 




'be 














0' 


be 







be 














ac 


Tk = 








ac 











ac 











ae 








ab 
















ab 





ah 



















ab 



Then, the local stiffness matrix Mk can be computed using the formula presented in 
Lemma 13.11 



4. Numerical Experiments. In this section, we shall report some numerical 
results for the weak Galerkin finite element method on a variety of testing problems, 
with different mesh and finite elements. To this end, let Uh — {uq, u;,} and u be 
the solution to the weak Galerkin equation (|2.2|) and the original equation (jl.ip . 
respectively. Define the error by Ch — Uh — QhU — {eo, Ch\ where QhU is the 
projection of u onto appropriately defined spaces. Let us introduce the following 
norms: 



semi-norm: 



Element-based norm : 



Edge/Face-based 1? norm : 



lleoll 




1/2 



Vde/J dx 



1/2 



Col dx 



where in the definition of ||ef,||, Kk stands for the size of the element K that takes F 
as an edge/face. We shall also compute the error in the following metrics 



llVdW?,- Vull 



1/2 



\IuV dx 



E 



1/2 



1^0 



u\ dx 



K 



leollc 



sup |eo(x)|. 

X e Ko 



Here the maximum norm ||eo||oo is computed over all Gaussian points, and all other 
integrals are calculated with a Gaussian quadrature rule that is of high order of 
accuracy so that the error from the numerical integration can be virtually ignored. 

4.1. Case 1: Model Problems with Various Boundary Conditions. First, 
we consider the Laplace equation with nonhomogeneous Dirichlet boundary condition: 



(4.1) 



u ^ g on dfl. 
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We introduce a discrete Dirichlet boundary data g^, which is either the usual nodal 
value interpolation, or the LP' projection of m = g on the boundary. Let T C and 
define 

Sg^X,h ={v ■■ v\k, e Pj{Ko) or QjiKo) for all K G 
v\FePi{F) or Qi{F) for all F e Fh, 
V ^ gh Fh'^ r}. 

When r = 9ri, we simply denote Sgi^^g^ji by Sgi^ji. The discrete Galerkin formulation 
for the nonhomogeneous Dirichlet boundary value problem can be written as: find 
Ufi G Sgi^^h such that for all u/i G 5"°, 

(AVdUfi.'^dVh) + (/3 ■ '^dUh,vo) + (7U0, vo) = if, vq). 

We would like to see how the weak Galerkin approximation might be affected when 
the boundary data u = g is approximated with different schemes (nodal interpolation 
verses L"^ projection). To this end, we use a two dimensional test problem with domain 
n — (0, 1) X (0, 1) and exact solution given by w = sin(27r2: + 7r/2) sin(27ry + 7r/2). 
A uniform triangular mesh and the element {Po{K), Po{F), RTo{K)) is used in the 
weak Galerkin discretization. The results are reported in Table 14.11 and Table 14.21 It 
can be seen that both approximations of the Dirichlet boundary data give optimal 
order of convergence for the weak Galerkin method, while the projection method 
yields a slightly smaller error in ||eo|| and ||ef,||. 

Next, we consider a mixed boundary condition: 

\u = g^ ouTd, 
1 (AVu) ■ n + au ^ g^ on Tr, 

where g^ is the Dirichlet boundary data, g^ is the Robin type boundary data, a > 0, 
and Tu n Tn = 0, U r^; = d^l. When a = 0, the Robin type boundary condition 
becomes the Neumann type boundary condition. 

For the mixed boundary condition, it is not hard to see that the weak formulation 
can be written as: find m/j G SgO such that for all u/j e Sqxo.h-, 

{AV dUh,y dVh) + {aub,Vb)rn + (/3 • V<ju/i,wo) + (7^0, vo) = (/,vo) + {9^,Vb)rR, 

where (•, ■)-ra denotes the inner-product on Fj^. We tested a two-dimensional 
problem with A to be an identity matrix and il = (0, 1)^ with a uniform triangular 
mesh. The exact solution is chosen to be m = sin(7ry)e~^. This function satisfies 

Vu ■ n -|- ii = 

on the boundary segment x — 1. We use the Dirichlet boundary condition on all other 
boundary segments. The element {Pq{K), Po{F), RTq{K)) is used in the discretiza- 
tion. For the Dirichlet boundary data, the LP projection is used to approximate the 
boundary data gj^ . The results are reported in Table 14.31 It shows optimal rates of 
convergence in all norms for the weak Galerkin approximation with mixed boundary 
conditions. 
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Table 4.1 

Case 1. Numerical results with Dirichlet data being approximated by the usual nodal point 
interpolation. 



h 


llVde/JI 


lleoll 


lleoll 


\\VdUh-yu\\ 


||mo - ^t|| 


lleolloo 


1/8 


7.14e-01 


2.16C-02 


4.05C-02 


l.Olc+O 


1.30C-01 


4.43C-02 


1/16 


3.56C-01 


5.61e-03 


l.Olc-02 


5.04G-01 


6.53e-02 


1.12e-02 


1/32 


1.78e-01 


1.41e-03 


2.53e-03 


2.51C-01 


3.27e-02 


2.86e-03 


1/64 


8.90e-02 


3.55e-04 


6.32e-04 


1.25C-01 


1.63e-02 


7.15e-04 


1/128 


4.45e-02 


8.88e-05 


1.57e-04 


6.29e-02 


8.18e-03 


1.79e-04 


r = 


1.0012 


1.9837 


2.0014 


1.0024 


0.9984 


1 QS7Q 


Table 4.2 

Case 1. Numerical results with Dirichlet data being approximated by 


projection. 


h 


llVde/JI 


lleoll 




\\VdUh-Vu\\ 


||mo - u\\ 


Ileolloo 


1/8 


7.10C-01 


1.75C-02 


3.08C-02 


l.Olc+0 


1.29e-01 


3.68C-02 


1/16 


3.55e-01 


4.59e-03 


7.69e-03 


5.04e-01 


6.52e-02 


9.54e-03 


1/32 


1.78e-01 


1.16e-03 


1.92e-03 


2.51e-01 


3.27e-02 


2.39e-03 


1/64 


8.90e-02 


2.90e-04 


4.81e-04 


1.25G-01 


1.63e-02 


6.01e-04 


1/128 


4.45C-02 


7.27e-05 


1.20C-04 


6.29e-02 


8.18e-03 


1.50e-04 


r = 


0.9993 


1.9808 


1.9999 


1.0015 


0.9968 


1.9861 



4.2. Case 2: A Model Problem with Degenerate Diffusion. We consider 
a test problem where the diffusive coefScient A is singular at some points of the 
domain. Note that in this case, the usual mixed finite element method may not 
be applicable due to the degeneracy of the coefficient. But the primary variable 
based formulations, including the weak Galerkin method, can still be employed for a 
numerical approximation. 

More precisely, we consider the following two-dimensional problem 

-V-{xyVu) = f infl, 
u — on dil, 

where = (0, 1)^. Notice that the diffusive coefhcient A = xy vanishes at the origin. 
We set the exact solution to be u = x{l — x)y{l — y). The configuration for the finite 
element partitions is the same as in test Case 1. We tested the weak Galerkin method 
on this problem, and the results are presented in Table 14.41 and Figure 14.11 

Since the diffusive coefficient A is not uniformly positive definite on $7, we have 
no anticipation that the weak Galerkin approximation has any optimal rate of conver- 
gence, though the exact solution is smooth. It should be pointed out that the usual 
Lax-Milgram theorem is not applicable to such problems in order to have a result 
on the solution existence and uniqueness. However, one can prove that the discrete 
problem always has a unique solution when Gaussian quadratures are used in the 
numerical integration. Interestingly, the numerical experiments show that the weak 
Galerkin method converges with a rate of approximately 0{h^'^) in || Vrfe;i||, 0(/i^-^^) 
in 1 1 Co 1 1 and ||e;,||. It is left for future research to explore a theoretical foundation of 
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Table 4.3 

Case 1. Numerical results for a test problem with mixed boundary conditions, where a Robin 
type boundary condition is imposed on part of the boundary. 



h 


l|Vde,|| 


Ikoll 


lle^ll 




||uo 


lleolloo 


1/8 


1.55C-01 


3.18C-03 


1.14e-02 


1.95C-01 


4.51C-02 


1.12C-02 


1/16 


7.87e-02 


8.20e-04 


2.90e-03 


9.82e-02 


2.25e-02 


3.18e-03 


1/32 


3.94e-02 


2.06e-04 


7.29e-04 


4.92e-02 


1.12e-02 


8.40e-04 


1/64 


1.97e-02 


5.17e-05 


1.82e-04 


2.46e-02 


5.64e-03 


2.15e-04 


1/128 


9.87e-03 


1.29e-05 


4.56e-05 


1.23C-02 


2.82C-03 


5.46C-05 


Oihn 
r — 


0.9958 


1.9876 


1.9926 


0.9971 


1 nnni 


1.9262 


Case 2. Numerical results for a 


Table 4.4 
test problem with degenerate 


diffusion A in 


the domain. 


h 


l|Vde/.|| 


Ikoll 


\\eb\\ 


||VdM/i - Vu| 


||uo 


lleolloo 


1/8 


5.61e-02 


3.32e-03 


6.60e-03 


5.75e-02 


5.48e-03 


1.27e-02 


1/16 


4.03e-02 


1.38e-03 


2.81e-03 


4.09e-02 


2.59e-03 


4.90e-03 


1/32 


2.95e-02 


5.68e-04 


1.16e-03 


2.96e-02 


1.23e-03 


2.21C-03 


1/64 


2.15e-02 


2.35e-04 


4.83e-04 


2.15e-02 


5.97e-04 


1.16e-03 


1/128 


1.55e-02 


9.93e-05 


2.02e-04 


1.55e-02 


2.91e-04 


5.99e-04 


om 

r = 


0.4614 


1.2687 


1.2594 


0.4697 


1.0579 


1.0912 



the observed convergence behavior. 

4.3. Case 3: A Model Problem on a Domain with Corner Singularity. 

We consider the Laplace equation on a two-dimensional domain for which the exact 
solution possesses a corner singularity. For simplicity, we take ft = (0,1)^ and let the 
exact solution be given by 

(4.2) uix,y)^x{l-x)y{l-y)r-^+\ 

where r — y/ x'^ + y^ and 7 G (0, 1] is a constant. Clearly, we have 

u e H^{n) n H^+-'-'{n) and u i 

where e is any small, but positive number. Again, a uniform triangular mesh and the 
element (Pq{K)., Po{F), RTq{K)) are used in the numerical discretization. Note that 
the weak Galerkin for this problem is exactly the same as the standard mixed finite 
element method. 

This model problem was numerically tested with 7 = 0.5 and 7 = 0.25. The 
convergence rates are reported in Table |475] and Table l476l Notice that ||Vde/i|| and 
||eo|| behaves in a way as predicted by theory (|2.4I) : i.e., they converge with rates given 
by 0{h'^) and 0{h^'^^), respectively. The result also shows that the approximation 
on the element edge/face has a rate of convergence O ( ft. ^ +'•'). 

4.4. Case 4: A Model Problem with Intersecting Interfaces. This test 
problem is taken from [TB] , which has also been tested by other researchers |T71 [20] . 
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Table 4.5 

Case 3. Convergence 'rates for a 'p'rohlcm w/fJi cornc/f si-n(p.dai-ii ij = O.o^. 



h 


W'^dehW 


l|eo|| 


I|e6|| 


\\VdUh - Vm|| 


\\uo-u\\ 


eo| oo 


1/8 


1.88e-01 


6.40e-03 


1.47e-02 


2.54e-01 


1.49C-02 


4.30C-02 


1/16 


1.36C-01 


2.20C-03 


5.28C-03 


1.84C-01 


7.66O-03 


3.01C-02 


1/32 


!).74r-()2 


7.62(--()4 


1.8()<--()3 


1.32o-()l 


3.89o-()3 


2.12r-()2 


1/64 


6.93C-02 


2.65C-04 


6.57C-04 


9.42C-02 


1.96C-03 


1.49C-02 


1/128 


4.92C-02 


9.33C-05 


2.32C-04 


6.69O-02 


9.88C-04 


1.05C-02 


0(/i'-) 
r = 


0.4852 


1.5251 


1.4992 


0.4827 


0.9805 


0.5066 



In two dimension, consider Q = (—1, 1)^ and the following problem 

-V • {AWu) = 0, 



where A = K1I2 in the first and third quadrants, and 7^212 in the second and forth 
quadrants. Here I2 is the 2x2 identity matrix and Ki , K2 are two positive numbers. 
Consider an exact solution which takes the foUowing form in polar coordinates: 



u{x,y) = r'^jjLie), 



where 7 G (0, 1] and 



(4.3) = { 



cos((7r/2 - ct)7) cos((6' - 7r/2 + p)7), 
C0s(p7) C0S{{9 — TT + cr)7), 
C0S((T7) C0S((^ — TT — p)7), 
cos((7r/2 - p)7) cos((6l - 37r/2 - a)^), 
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if < 6* < 7r/2, 
if 7r/2 < 6* < TT, 
if TT < (9 < 37r/2, 
if 37r/2 <e<2-K. 



Table 4.6 

Case 3. Convergence rates for a problem with corner singularity (7 = 0.25^. 



h 




Ikoll 




\\VdUh-Vu\\ 


||iio 


lleolloo 


1/8 


4.93e-01 


1.69e-02 


3.58e-02 


6.65e-01 


2.56e-02 


1.25C-01 


1/16 


4.18e-01 


7.07e-03 


1.52e-02 


5.66e-01 


1.31e-02 


1.05e-01 


1/32 


3.53e-01 


2.94e-03 


6.39e-03 


4.79e-01 


6.72e-03 


8.85e-02 


1/64 


2.98e-01 


1.22e-03 


2.68e-03 


4.04e-01 


3.42e-03 


7.44e-02 


1/128 


2.51e-01 


5.14C-04 


1.12G-03 


3.40e-01 


1.73e-03 


6.25e-02 


O(ft') 
r = 


0.2437 


1.2613 


1.2489 


0.2417 


0.9717 


0.2505 



The parameters 7, p, a satisfy the following nonhnear relations 

R := K1/K2 = - tan((7r/2 - (7)7) cot(p7), 
1/i? = — tan(p7) cot((T/9), 
(4.4) - tan(p7) cot((7r/2 - p)7), 

max{0, 7r7 — tt} < 27/3 < min{7r7, tt}, 
max{0, TT — 7r7} < < min{7r, 2ti — 717}. 

The solution u{r, 9) is known to be in H'^^^'^'^ {^) for any e > 0, and has a singularity 
near the origin (0, 0). 

One choice for the coefficients is to take 7 = 0.1, R « 161.4476387975881, 
p « 7r/4, a K, —14.92256510455152. We numerically solve this problem by using the 
weak Galerkin method with element {Pq{K)^ Po{F), RTq{K)) on triangular meshes. 
It turns out that uniform triangular meshes are not good enough to handle the sin- 
gularity in this problem. Indeed, we use a locally refined initial mesh, as shown in 
Figure which consists of 268 triangles. This mesh is then uniformly refined, by di- 
viding each triangle into 4 subtriangles, to get a sequence of nested meshes. Although 
this can not be compared with an adaptive mesh refinement process, it does improve 
the accuracy of the numerical approximation, as shown in our numerical results re- 
ported in Table l477l Since the mesh is not quasi- uniform, we do not expect that the 
theoretic error estimation (|2.4p apply for this problem. An interesting observation 
of Table H771 is that, the norm ||uo — u|| appears to converge in a much faster rate 
than ||eo|| — ||uo — Qo^illi while the opposite has usually been observed for other test 
cases. We believe that this is due to the use of a locally refined initial mesh in our 
testing process. When the actual value of ||uo — u\\ reduces to the same level as the 
value of Ijeoll, its convergence rate slows down to the same as ||eo||. Readers are also 
encouraged to derive their own conclusions from these numerical experiments. 

We also observe that, when the initial mesh gets more refined near the origin, 
the convergence rates increase slightly. In Table 14. 8[ this trend is clearly shown. 
For each initial mesh, it is refined four times to get five levels of nested meshes. The 
convergence rates are computed based on these five nested meshes. The initial meshes 
are generated by refining only those triangles near the origin. Two examples of initial 
meshes are shown in Figure 

4.5. Case 5: An Anisotropic Problem. Consider a two dimensional anisotropic 
problem defined in the square domain = (0, 1)^ as follows 



-V • {AVu) = /, 
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Fig. 4.2. Case 4- The initial triangular mesh for the intersecting interface problem, with 268 
(left) and 300 (right) triangles. 




Table 4.7 

Case 4- Convergence rate for the intersecting interface problem with an initial mesh containing 
268 triangles. 



level 




Ikoll 




\\VdUh-yu\\ 


||wo - "11 


||eo||oo 





1.07e-01 


3.97e-03 


9.95e-03 


1.47e-01 


2.60e-02 


1.97C-02 


1 


9.76e-02 


2.92e-03 


6.44e-03 


1.26C-01 


1.33e-02 


1.94e-02 


2 


9.30C-02 


2.51e-03 


5.11e-03 


1.16e-01 


7.01e-03 


1.91C-02 


3 


9.12C-02 


2.21e-03 


4.44e-03 


l.llc-01 


3.95e-03 


1.88e-02 


4 


8.98C-02 


1.95e-03 


3.91e-03 


1.07e-01 


2.55e-03 


1.84e-02 


r = 


0.0604 


0.2446 


0.3229 


0.1084 


0.8461 


0.0239 



where the diffusive coefScient is given by 



A 



P 6 

1 



for fc ^ 0. 



We chose a function / and a Dirichlet boundary condition so that the exact solution 
is given by u{x,y) = sin(27ra;) sin(2fc7r?/). In applying the weak Galerkin method, we 
use an anisotropic triangular mesh that was constructed by first dividing the domain 
into kn X n sub-rectangles, and then splitting each rectangle into two triangles by 
connecting a diagonal line. The characteristic mesh size is h = 1/n. We tested two 
cases with fc = 3 and k — 9. The results are reported in Tables 14.91 and 14.101 The 
tables show optimal rates of convergence for the weak Galerkin approximation in 
various metrics. The numerical experiment indicates that the weak Galerkin method 
can handle anisotropic problems and meshes without any trouble. 

4.6. Case 6: A Three-Dimensional Model Problem. The final test prob- 
lem is a three dimensional Laplace equation defined on i7 = (0, 1)'^, with a Dirichlet 
boundary condition and an exact solution given by m = sin(27ra;) sin(27r2/) sin(27r2;). 
The purpose of this test problem is to examine the convergence rate of the cubic 
(QoiK), QoiF), RTq{K)) element. The results are reported in Table SUTl 

In addition to the optimal rates of convergence as shown in Table HTTTl on can also 
see a superconvergence for ||Vde?i||. The same result is anticipated for 2D rectangular 
elements. It is left to interested readers for a further investigation, especially for 
model problems with variable coefficients. 
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Table 4.8 

Case 4- Convergence rate for the intersecting interface problem with different initial meshes, 
where the first column indicates the total number of triangles in the initial mesh. 



41 

it 

triangles 




Convergence rates 0{h^), r 
''() ' V,/(//, - V(/ 


"0 - " 






268 


0.0604 


0.2446 


0.3229 


0.1084 


0.8461 


0.0239 




300 


U.U ( OU 


0.2623 


0.3489 


0.1206 


U.soyy 


U.Uo ( 6 




332 


0.0888 


0.2818 


0.3772 


0.1329 


0.8912 


0.0487 




364 


0.1020 


0.3031 


0.4079 


0.1454 


0.9099 


0.0586 




396 


0.1148 


0.3266 


0.4411 


0.1581 


0.9260 


0.0673 




428 


0.1273 


0.3522 


0.4766 


0.1711 


0.9396 


0.0749 




460 


0.1396 


0.3802 


0.5145 


0.1843 


0.9509 


0.0817 




492 


0.1519 


0.4105 


0.5548 


0.1978 


0.9602 


0.0878 




524 


0.1641 


0.4432 


0.5972 


0.2117 


0.9678 


0.0932 








Case 5. 


Table 4.9 

Convergence rate for the anisotropic probl 


sm with k = 


3. 




h 


l|Vde;.|| 


lleoll 


l|e.|| 


IIVdMh - Vm|| 


\\un — u\\ 

1 1 --^u 1 1 


||eo||oo 


1/8 


1.48e+0 


1.95C-02 


4.61e-02 


2.70e+0 


1.29C-01 


4.13C-02 


1/16 


7.39e-01 


5.11C-03 


1.16e-02 


1.35e+0 


6.53C-02 


1.06C-02 


1/32 


3.69e-01 


1.29e-03 


2.92e-03 


6.80e-01 


3.27e-02 


2.67e-03 


1/64 


1.84e-01 


3.24e-04 


7.33e-04 


3.40e-01 


1.63e-02 


6.68e-04 


1/128 


9.23e-02 


8.12e-05 


1.83e-04 


1.70e-01 


8.18e-03 


1.66e-04 


r = 




1.0010 


1.9793 


1.9942 


0.9972 


0.9975 


1.9906 






Case 5. 


Table 4.10 
Convergence rate for the anisotropic probl 


;m with k = 


9. 




h 




V,/(';m 


ikoll 






"0 - " 


'■Ill 


1/4 


7.98e+0 


6.80C-02 


2.93C-01 


1.58e+l 


2.52C-01 


1.49e-01 


1/8 


3.89e+0 


2.07e-02 


7.44e-02 


8.18e+0 


1.30e-01 


4.22e-02 


1/16 


1.91e+0 


5.43e-03 


1.88e-02 


4.12e+0 


6.53e-02 


1.09e-02 


1/32 


9.54C-01 


1.37e-03 


4.72C-03 


2.06e+0 


3.27C-02 


2.74C-03 


1/64 


4.76e-01 


3.44e-04 


1.18e-03 


1.03e+0 


1.63e-02 


6.84e-04 


)■ = 




I.OIGI 


1.9160 


1.9897 




0.988.3 


1.9492 


Table 4.11 

Case 6. Convergence rate for a 3D model problem with smooth solution. 


h 




Vde/j|| 


l|eo|| 


\\eb\\ 


llVdUfe - Vu|| 


||uo ~ u\\ 


eo cx) 


1/8 


1.85C-01 


1.62C-02 


4.27C-02 


1.22C+00 


1.34C-01 


3.63C-02 


1/12 


8 


.53e-02 


7.69e-03 


1.94e-02 


8.19e-01 


9.14e-02 


1.96e-02 


1/16 


4.86e-02 


4.42e-03 


l.lOe-02 


6.15e-01 


6.89e-02 


1.18e-02 


1/20 


3.13e-02 


2.85e-03 


7.07e-03 


4.92e-01 


5.52e-02 


7.78e-03 


o{hn 

r = 




1.9389 


1.8984 


1.9618 


0.9914 


0.9737 


1.6779 
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